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ABSTRACT 

An enhanced rate of stellar tidal disruption events (TDEs) may be an important 
characteristic of supermassive black hole (SMBH) binaries at close separations. Here 
we study the evolution of the distribution of stars around a SMBH binary due to 
the eccentric Kozai-Lidov (EKL) mechanism, including octupole effects and apsidal 
precession caused by the stellar mass distribution and general relativity. We identify a 
region around one of the SMBHs in the binary where the EKL mechanism drives stars 
to high eccentricities, which ultimately causes the stars to either scatter off the second 
SMBH or get disrupted. For SMBH masses 1O 7 M 0 and 1O 8 M 0 , the TDE rate can 
reach ~ 10 _2 /yr and deplete a region of the stellar cusp around the secondary SMBH 
in ~ 0.5 Myr. As a result, the final geometry of the stellar distribution between 0.01 
and 0.1 pc around the secondary SMBH is a torus. These effects may be even more 
prominent in nuclear stellar clusters hosting a supermassive and an intermediate mass 
black hole. 
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1 INTRODUCTION 

Supermassive black holes (SMBHs) are ubiquitous at the 


centers of galaxies (Kormendy & Ho 20131. Stars passing 


close to the SMBH can be tidally disrupted, and the fall 
back of the stellar debris produces a strong electromagnetic 
tidal disruption flare (e.g., |Gezari 2012). More than a dozen 


tidal disruption event (TDE) candidates have been observed 


until present (e.g., 

Bade et al.|1996 

Gezari et al. 2003 

2006 

to 

o 

o 

OO 

2009; (van Ye 

zen et al.|2011 Gezari et al.[2012; 1 

oloien 


et al. 20141, including two candidates with relativistic jets 
( Levam et al.|20lT Bloom et al.||20lT Zauderer et al.|2011 


Cenko et al.|2012 i. TDEs can provide valuable information 
on dormant SMBHs, which are otherwise difficult to detect. 

The rate of the TDEs provide information about the 
SMBH and the stellar distribution in the center of galaxies 
( |Stone fe Metzger]|2014| . The rate of TDEs is highly un¬ 
certain observationally due to the small sample size. It is 
estimated to be in the range of 10 -5 — 10 -4 per galaxy 
per year by |Donley et ah] (|2002[) ; |Gezari et al7| ( |2008[ 1; 
Maksym (2012l7|v an Velzen fc Farrar (20141. This roughly 


agrees with the theoretical estimates, discussed by Frank 

& Rees (1976); Lightman & Shapiro) (jl977|); |Cohn & Kul- 

srud 

(1978 

; Magorrian & Tremaine (19991; Wang & Mer- 

ritt l 

o 

o 

CM 

Brockamp et al. (2011); Stone & Metzger (2014). 


However, the TDE rate may be enhanced due to the pres¬ 
ence of a non-axisymmetric gravitational potential around 


the SMBH (M erritt fc Poon|200 4;), or due to a massive per- 
turber ( Perets et al.||2007 |. In addition, the TDE rate may 
be higher in galaxies with more than one SMBH (|Ivanov| 


et al.|[2005 Chen et ah|2009[ 2011| Wegg fc BodepoTT I, or 

when the SMBH binary (SMBHB) recoils due to the emis¬ 
sion of gravitational waves ( |Stone fc Loeb|[2011| |Li et alT 
2012 Stone fc Loeb|2012|. Some TDEs may not appear as 


flares and therefore be missed in observations! Guillochon & 
Ramirez-Ruiz|2015 I. 


The interaction between a SMBHB and an ambient star 
cluster has been discussed in the literature using numerical 


scattering experiments by Sesana et al. (2011I and using di¬ 


rect N-body simulations by Iwasawa et al. (20111; Gualan 


dris & Merritt (20121; Meiron & Laor (20131; Wang et al. 


(20141. In particular, it has been shown that the star clus¬ 


ter may either increase or decrease the eccentricity of the 
SMBHB depending on the fraction of counter-rotating to 
co-rotating stars. The SMBHB ejects a population of stars 
from the cluster in an anisotropic manner, and the SMBHB 
produces a deficit in the number density of stars, a dip in 
the velocity dispersion in the inner regions, and an inner 
counter-rotating and an outer co-rotating torus of stars with 
respect to the binary. 

In this paper, we focus on the distribution of stars or¬ 
biting close to one member of the SMBHB and perturbed 
by the other SMBH through hierarchical three-body inter- 
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actions. We examine the effect of these hierarchical three 
body interactions. Specifically, the outer SMBH perturbs 
the stellar population around the inneiQ SMBH, and leads 
to long-term variations in the eccentricities and inclinations 
of the stellar orbits while keeping the semimajor axes of their 
orbits fixed. In particular, when the orbit of the SMBH sec¬ 
ondary is circular and if the mutual inclination between the 
orbits of the SMBHB and a star is over 40°, the stellar eccen¬ 
tricity and inclination undergo periodic oscillations, known 
as the quadrupole Kozai-Lidov mechanism (Kozai|1962| Li 


dov|1962|. This is caused by the long-term (orbit-averaged) 


Newtonian gravitational effect expanded in multipoles to the 
quadrupole order, i.e. second order in the semimajor axis ra¬ 
tio of the stellar and the outer SMBH’s orbit. More generally, 
it has been found that when the outer orbit is eccentric, the 
analogous octupole eccentric Kozai-Lidov mechanism (EKL, 
third order in semi major axis ratio) causes the eccentricity 
to be excited very close to unity and the inner orbit to flip 
relative to the invariable plane from prograde to retrograde 


or vice versa (e.g., 

Ford et al. 

2000 

Naoz et al. 

2011 fEati 

et af. 2011. Lithwick & Naoz 

2011 

|Naoz et al. 

013a b ; Li 

et al. 014a|b). The TDE rate has been discussed in the lit- 


erature for stars orbiting an SMBHB, where the quadrupole 
Kozai-Lidov mechanism can enhance the tidal disruption 
rate (llvanov et ~ah 2005 Chen et al. 2009 2011| Wegg & 


Bode 2011). For the Galactic Center, the Kozai-Lidov mech¬ 


anism driven by the stellar disk has also been discussed and 
the additional effects of Newtonian apsidal precession were 
shown to play a significant role ( Chang|200 9). In light of re¬ 
cent developments in the understanding of hierarchical three 
body interactions we revisit this problem. Since the stellar 
eccentricity can be increased to a value much closer to unity 
by eccentric perturbers, we expect the EKL mechanism to 
enhance TDE rates with respect to the circular case. We 
therefore seek to re-evaluate the total number of stars vul¬ 
nerable to TDE due to EKL. 

It is well known that apsidal precession quenches the 


EKL mechanism (e.g., Ford et al. 2000 Blaes et al. 2002 
Naoz et al.|013b'|. In galactic nuclei this may be due to the 


Newtonian (NT) gravitational effect of the spherical stellar 
cusp or general relativistic (GR) precession, provided that 
the corresponding precession timescale is much shorter than 
the Kozai timescale (Chang||~2009). Furthermore, the EKL 
mechanism may be quenched if the eccentricity of the star 
is changed by the stellar cluster due to scalar resonant re¬ 
laxation, or if the orbital plane is reoriented by the stellar 


cluster due to vector resonant relaxation (Rauch & Tremaine 
1996 Kocsis fe Tremaine|2011 20141 or Lense-Thirring pre- 


cession ( |Merritt et al. ||2010| |Merritt fc Vasiliev||2012| ). We 
find that NT precession and GR precession may have a large 
effect on the EKL mechanism, but tidal effects, scalar and 
vector resonant relaxation, and Lense-Thirring precession 
are typically less important. The timescale on which the 
EKL mechanism operates increases if the outer SMBH mass 
is reduced. Thus, GR precession may dominate over and 
quench the EKL mechanism most efficiently if the outer 
SMBH is less massive than the inner SMBH (see figure 2 


1 We consider stars that initially orbit the “inner” SMBH and 
whose orbits are perturbed by the “outer” SMBH regardless of 
which SMBH is more massive. 


Naoz & Silk 2014). Similarly, we find that NT preces¬ 


sion also suppresses the EKL mechanism most efficiently 
when the outer SMBH is less massive. Tidal disruption is 
expected in the opposite regime when the EKL mechanism 
is very prominent, i.e. when the outer SMBH is more mas¬ 
sive than the inner SMBH. We identify the outcome of the 
EKL mechanism as a function of SMBHB parameters and 
quantify the TDE rate. 

Our discussion is organized as follows. In §2, we describe 
the adopted methods. In §3, we characterize the parameter 
space to identify where the EKL mechanism is important. 
Then, we calculate the tidal disruption rate and discuss the 
final stellar distribution due to the EKL mechanism with 
an illustrative example in §4, and for stars surrounding an 
intermediate-mass black hole in §5. Finally, we summarize 
our main results in §6. 


2 METHOD 


We study the tidal disruption of stars due to the EKL mech¬ 
anism in galaxies that host a SMBHB. The three-body sys¬ 
tem consists of an “inner binary” comprised of the SMBH 
and a star, and an “outer binary” comprised of the outer 
SMBH and the center mass of the inner binary, as shown in 
Figurejl] We denote the masses of the objects by mo (inner 
SMBH), mi (star), and m 2 (outer SMBH), and for orbital 
parameters we use subscript 1 and 2 for the inner and outer 
binary, respectively. In order for the EKL mechanism to op¬ 
erate, we require the triple system to be in a hierarchical 
configuration: the inner binary on a much tighter orbit than 


the third object, such that (e.g., Lithwick & Naoz 2011 Katz 
et al.|2011| >, 


a 1 e2 


< 0.1, 


(1) 


where a and e are respectively the semimajor axis and ec¬ 
centricity. 



Figure 1. The system configuration, ‘c.m.’ denotes the center of 
mass of the inner binary, which contains the star (with mass mi) 
and SMBH (with mass mo). The other SMBH (with mass m 2 ) is 
on an outer orbit. 
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2.1 Comparison of timescales 

We examine the range of orbital parameters in oder to iden¬ 
tify the regions in which the EKL mechanism may operate. 
The relevant processes’ timescales can be expressed as: 


27ra 3 (l - e|) 3/2 x/(m 0 + mi)(l - e?) , 

tK = - ^ I/O - \ Z ) 


\[Gaf 2 mi 


toct — tK 

e 


27ra 3/ V(l - e{) 
3G 3 / 2 (m 0 + mi) 3 / 2 


ton 


2na 5 2 /2 c 2 (l - e|) 

GR ' J 3G 3 / 2 (mo + mi + ms) 3 / 2 
, 16 a\c 2 {\-e\fl 2 (mof 2 

tGR,int — ~ 


t NT — 27T 


9 y/aiei\fl — e\G z ! 2 rn\rn2 

p /Gmo/af 


tRR,s = 


AtOjJ 


7rmoei 

ml 


- J dip M» (r) cos ip'j 


t3 2 Q. 2 M*(r)mi 


2nf vrr 

tRR,v = - 71 - m 0 


~ / / f-u , 

“ ^/M,(r)mi 

_3 


Cei = 0.34 


t.LT = 


tGW = 


G 2 pm i In A 
afc 3 (l — e 2 ) 3 / 2 


2G 2 r 


ai 5 


4 64 G 3 mom2(mo + m 2 ) 


( 3 ) 

( 4 ) 

( 5 ) 

( 6 ) 

( 7 ) 

( 8 ) 
( 9 ) 

( 10 ) 

( 11 ) 

( 12 ) 


Here tx is the quadrupole (C(ai/o 2 ) 2 ) Kozai timescale. Fol¬ 


lowing Naoz et al. (013b), t oc t is the octupole (0(ai/a2) 3 ) 
Kozai timescale. tGRi and tGR 2 are the timescales of the 
first order post Newtonian general relativistic (GR) preces¬ 
sion at the quadrupole order ( 0 (ai/a 2 ) 2 ) on the inner and 
outer orbit, and tGR,int is the timescale associated with the 
first post-Newtonian order GR interaction between the inner 


and the outer orbit. Following Kocsis & Tremaine (20111, 


t nt is the timescale of the Newtonian precession caused by 
the stellar potential, and fxK, s and txx.u are the timescales 
of the scalar and vector resonant relaxation. t re i is the two 
body relaxation timscale. t lt is the Lense-Thirring preces¬ 
sion timescale, and tGW is the timescale of the orbital decay 
of the binary SMBHB due to gravitational wave radiation. 
For the resonant relaxation timescales, M* (r) is the mass of 
the stars interior to r, u is the net rate of precession due 
to GR and NT, /3 S is estimated to be 1.05 ± 0.02 by |Eilon] 
et al. (2009), H is the orbital frequency of the star, and f vrr 


is estimated to be 1.2 by Kocsis & Tremaine (2014 1 . For 
the Lense-Thirring timescale, sGml/c is the spin angular 


Naoz 


momentum of the inner SMBH (see references in e.g. 
et al.|013b Kocsis &; Tremaine|2011 Peters|1964 1. We define 
some of these effects in more detail in § |2.2| below. 


The EKL mechanism operates if the following criteria 
are satisfied: 

(i) The three-body configuration satisfies the hierarchical 
condition (e < 0.1, see equation 0 ) 

(ii) The stars stay in the Hill sphere of the inner SMBH 


in order for them to remain bound to it, i.e. ai(l + ei) < 
02(1 - e 2 )(mo/3m 2 ) 1 7 3 . 

(iii) The quadrupole (G(ai/a 2 ) 2 ) Kozai timescale, tx, 
needs to be shorter than the timescales of the other mecha¬ 
nisms that modify the orbital elements, otherwise the EKL 
mechanism is suppressed. The competing mechanisms in¬ 
clude Newtonian precession (NT), general relativistic pre¬ 
cession (GR), scalar resonant relaxation, vector resonant 
relaxation, two-body relaxation, Lense-Thirring precession, 
and the gravitational radiation. 

Note that the secular approximation fails when the per¬ 
turbation from the outer SMBH is too strong or when the 


eccentricity reaches values very close to unity (e.g., Antonini 
fe Perets||2012 |Katz fe Dong||2012 |Antognini et al.||2014| 


Antonini et al. 2014 Bode & Wcgg 20141. This means that 
there are some systems that are poorly described by our ap¬ 
proximation. However, we expect that those systems reach 
even higher eccentricities than the one predicted by the oc¬ 
tupole approximation (e.g., Antognini et al. 2013), and thus 
our overall qualitative conclusions may hold even for those 
systems, but the quantitative rate values possibly underes¬ 
timate the true rates. 

To calculate the Newtonian timescale, the resonant re¬ 
laxation timescales, and the two body relaxation timescale, 
we adopt the spherically symmetric model for the stellar 


density discussed in O’Leary & Loeb (20091. Specifically, 


the stellar density distribution is a power law of semimajor 
axis and the normalization is fixed by the M — a relation, 


3 —a m 0 (GMotmo/Mo) 1 - 27 * 


\ 


-3+ex 


(13) 


where k = 4, Mo = 1.3 x 10 8 Mq, <to = 200km/s (Tremaine 
et al.|2002| , and we set a = 1.75. 

Figure p] shows the timescales for the case of a 1 M@ star 
orbiting a 10 7 M 0 SMBH. The separation of the SMBHB is 
set to 0.3 pc. The upper panel corresponds to m 2 = 10 G M@, 
and the lower panel corresponds to m 2 = 10 9 Mg. For the 
Lense-Thirring timescale, s is set to unity. The eccentric¬ 
ity of the star-SMBH system, ei, is assumed to be 2/3 
and e 2 is assumed to be 0.7. The EKL-dominated region 
is larger for higher e 2 with fixed ai and 02 . Figure [2] shows 
that the EKL mechanism is suppressed for a 1O 7 -1O G M 0 bi¬ 
nary at all radii, but it may operate at least in a restricted 
range for a 10' 10 9 M 0 binary. Note that although the oc¬ 
tupole timescale t oc t is longer than some of the other sec¬ 
ular timescales, our simulations show that the eccentricity 
can nevertheless reach high values provided that tx is the 
shortest timescale and t oc t is at most moderately larger than 
the other timescales. Since t oc t = tx/e and 1/e ~ 10 — 30, 
toct is only moderately larger than the other timescales in 
most of the relevant phase space when tx is the shortest 
timescale. Thus, in the following, we identify the regions 
where the eccentricity may be excited using conditions f/T) 


(iii) above irrespective of t oc t. Typically, the conditions on 


the quadrupole Kozai timescale (tx < tGR and tx < /jvt) 
set the lower limit for ai for a fixed 02 , and the hierarchical 
configuration e < 0.1 and the Hill sphere limit set the upper 
limit on ai. 

Next, we examine the <21 — 02 parameter space to iden¬ 
tify the parameters where EKL dominates. We plot two 
examples in Figure [ 3 ] mo = 10' M 0 , mi = 1M 0 , m 2 = 
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Figure 2. The different timescales as a function of the semi 
major axis of the stars (ai), where ei = 2/3, m o = 10 7 Mq, Q 2 = 
0.3 pc, mi = 1M 0 , e 2 = 0.7. In the upper panel, m 2 = 10 6 Mq, 
and in the lower panel, m 2 = 10 9 Mq. In the grey region, e > 0.1, 
the hierarchical approximation is violated. The EKL mechanism 
does not operate in the grey region and wherever is not 

the shortest timescale. The quadrupole Kozai timescale is shorter 
than the other timescales for the semimajor axis range indicated 
by the light green arrow. 

10 6 Mq, e 2 = 0.7 in the upper panel, and mo = 10' M@, 
mi = IMq, m 2 = 10 9 Mq, e 2 = 0.7 in the lower panel. 
The EKL-dominated region is bigger for larger e 2 - To test 
the dependence on ei, we include two ei values: ei = 0.001 
(solid lines) and ei =2/3 (dashed lines), where ei = 2/3 
corresponds to the mean value of ei in a thermal distribu¬ 
tion. The parameter space is independent of the mass of the 
star as long as mi <C mo. The EKL-dominated region is 
bounded by tx = tGR (blue line) and fx = tjvT (red line) 
from above and by the Hill sphere limit (grey line) and the 
hierarchical condition (black line) from below. In the upper 
panel, there is no region where the EKL mechanism domi¬ 
nates. In the lower panel, the region where EKL dominates 
is shaded with horizontal dashed lines for ei = 2/3 and it is 
shaded with vertical solid lines for ei = 0.001. 

We calculate the number of stars affected by the EKL 
mechanism for the particular stellar density distribution 
around the inner SMBH (equation m- In Figure [f] we 
consider the parameter space of different mo, m 2 , < 22 , e 2 and 
show the number of stars in the range of <21 where all cri¬ 
teria are satisfied for the EKL mechanism to operate. Each 
panel shows the parameter plane of mo and m 2 (assuming 
mi -C mo), 02 is varied in different columns of panels from 
0.1 to 10 pc, and e 2 is varied in the different rows from 
0.1 to 0.7. We set the stellar eccentricity to ei = 2/3 in all 
panels, the mean eccentricity for an isotropic thermal dis¬ 
tribution. There is no systematic change in the number of 
stars affected by the EKL versus e\. When ei = 0.001, the 




Figure 3. The ai — <22 parameter space, mo = 10 7 Mq, mi = 
IMq, e 2 = 0.7. In the upper panel, m 2 = 10 6 Mq, and in the 
lower panel, m 2 = 10 9 Mg. The solid blue and red lines represent 
ei = 0.001 and the dashed blue and red lines represent ei = 2/3. 
Above the red or blue lines, the EKL mechanism is suppressed by 
the GR or the Newtonian precession. Below the black line or the 
grey lines, the hierarchical configuration or the Hill sphere limit 
is violated. The EKL mechanism is suppressed everywhere in the 
upper panel, and the EKL mechanism dominates in the shaded 
regions in the lower panel. 

numbers typically increase to roughly twice the numbers of 
ei = 2/3, since the maximum <21 allowed due to the Hill 
sphere criterion becomes larger. When ei = 0.999, the pa¬ 
rameter region where stars can be affected in the mo — m 2 
plane increases, since the Newtonian precession timescale 
increases, while the changes in the numbers depend on the 
specific mo — m 2 configurations. In regions where the EKL 
mechanism is important, approximately 10 5-6 stars are af¬ 
fected. Thus, the EKL mechanism may significantly con¬ 
tribute to the tidal disruption events. Note that the EKL 
mechanism is more likely to be suppressed for stars orbiting 
around the more massive SMBH. However, for parameters 
where the EKL mechanism is not suppressed everywhere 
around the more massive inner SMBH, the total number of 
stars affected by EKL may be higher for stars orbiting the 
more massive SMBH than for those orbiting the less massive 
SMBH. 

2.2 Equations of motion 

As shown in the previous section, GR and NT precessions 
represent important limitations for the EKL mechanism. 
In this section, we review the equations of motion which 
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Figure 4. The number of stars (N) influenced by the EKL mechanism, assuming a stellar density distribution in equation < |13| ), and 
that the stellar mass is negligible and ei = 2/3 (the mean eccentricity assuming a thermal distribution). We determine the range of 
stellar semimajor axis ai where the EKL mechanism operates for a fixed set of SMBH masses, mo, m 2 , and outer orbital parameters, 
e 2 and 02 . Plotting the corresponding number of stars as a function of mo and m 2 for an array of e 2 and a 2 as shown, captures a large 
parameter space. The EKL mechanism affects a large number of stars over a wide range of SMBH binary parameters when 012 < 3 pc. 


govern the long-term evolution of stars due to the EKL 
mechanism, GR and NT precessions, and tidal effects 
adopted from Naoz et al. (013a|b I and Tremaine (20051. 


We use the Delaunay’s elements, which provide a con¬ 
venient dynamical description of hierarchical three-body 
systems. The coordinates are the mean anomalies, h 
and I 2 , the arguments of periastron, gi and g 2 , and the 
longitude of nodes, hi and /i 2 - Their conjugate momenta are 


Li = m ° mi yj G(mo + TOi)iii (14) 

mo + mi 

T m 2 (m 0 + mi) -■-■-r— 

L 2 = -V G (mo + mi + m 2 )a 2 

m 0 + mi + m 2 

Gi = LisJl-el,G 2 = L 2 \]l-el (15) 

Hi=Gicosii, H 2 =G 2 Cosi 2 , (16) 


where i denotes the inclination relative to the total angular 
momentum of the three-body system and G without sub¬ 
script is the gravitational constant. To leading order, the 
two binaries follow independent Keplerian orbits where lj 
are rapidly varying and Lj, Gj, Hj, gj, and hj are con¬ 
served for j £ {1,2}. These quantities are slowly varying 
over longer timescales due to the superposition of the per¬ 
turbations: the EKL mechanism, GR and NT precessions, 
and tidal effects, discussed next. 


2.2.1 Eccentric Kozai-Lidov Mechanism 


The equations of motion for the EKL mechanism may be 
derived using the double averaged Hamiltonian (i.e. aver¬ 
aged over the rapidly varying h and I2 elements). We go 


beyon d the analyses of Chen et al. (20111 and Wegg & Bode 
(20111, who considered only the quadrupole ( 0 (ai/a 2 ) 2 ) 
Kozai-Lidov mechanism, where the z-component of angu¬ 
lar momentum is constant. This assumption does not hold 
when the orbit of the SMBHB is eccentric, and one needs 


to include the octupole order terms (C?(ai/a 2 ) 3 ) (e.g. Naoz 
et al. |013a I. The Hamiltonian can be decomposed as 

d = G 2 {(2 + 3e?)(3cos 2 i tot - 1) 


H 


K ozai,quad 


+ 15e? sin” itot cos 2gi} 


(17) 


15 2 

LL k ozai,oct = — t M eiC 2 { A cos (j) + 10 cos it ot sin i to t 


x (1 — el) sin (/i sin < 72 }, 


(18) 


where 


mo — m\ 

6m — -;-€ 

m 0 + mi 

„ _ G 2 (mo + mi)' _ m^ Lj 

16 (mo + mi + m 2) 3 (momi) 3 L\G 3 

A = 4 + 3ei — ^ B sin 2 i tot 

B = 2 + 5e 1 — 7ei cos 2^i 
cos <j> = — cos gi cos <72 — cos itot sin gi sin g 2 


(19) 

( 20 ) 

( 21 ) 

( 22 ) 

(23) 
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The equations of motion for the EKL mechanism are given 
by Hamilton’s equations (eqn (A26-35) in Naoz et al.|013a). 


evolution with and without tides. The effect of tides is neg¬ 
ligible mainly because the orbital period is long and Q is 
low. 


2.2.2 GR effects 


Next, we consider the leading order (first Post-Newtonian, 
1PN) effects of GR. We follow Naoz et al. (013bI, who de¬ 
rived the double averaged 1PN Hamiltonian to the octupole 
(0(ai/a,2) 3 ) order. The Hamiltonian consists of four terms: 
Ha i, Ha 2 , H ai a 2 , Hint (Naozet al. 013b). Here H ai a 2 does 
not contribute to the dynamical evolution, and the long¬ 
term effect of Hint is typically negligible, as its timescale is 
longer than that of the Kozai timescale and the GR preces¬ 
sion of the inner and outer orbit as long as the star stays 
within the Hill sphere of the inner SMBH. Thus, we only 
consider the effects of H ai and Ha 2 which cause the GR 
precession of the arguments of periapsides, 


dg i I _ 3G 3/2 (m 0 + mi) 3/2 

dt\iPN, ai o 3/2 c 2 (1 — e 2 ) 

dg 2 \ _ 3 G 3/2 (m 0 + mi + rri2) 3/2 

dt \iPN,a 2 a 2 /2 c 2 (l-e|) 


(24) 

(25) 


Given that we neglect Hint, and higher order Post- 
Newtonian corrections such as Lense-Thirring precession 
and gravitational radiation, the other conserved quantities, 
Lj, Gj, Hj, hj, are not effected for j £ {1,2}. 


2.2.3 NT precession 


The Newtonian potential of a spherical stellar cusp causes 
apsidal precession at the rate (Tremaine|20051: 


gi,NT — 


(1 - el) 1 ' 2 




(Gmo/af) 1 / 2 aiei dr 


cos if, 


(26) 


where <f>* is the stellar potential, r is the distance to the cen¬ 
tral SMBH and ip is the true anomaly of the inner orbit. The 
averaged precession rate of gi due to Newtonian precession 
is expressed below: 


9i,nt 


Q Gmo/al) 1 / 2 

nmod 


rTT 

Jo 


dip M* (r) cos ip, (27) 


where A/*(r) is the mass of the stellar system interior to r 
and r = r(ip) = oi (1 — e 2 )/(l + e cos ip) from Kepler’s equa¬ 
tion. Explicit analytic expressions for the apsidal precession 
rate are given in Appendix A of Kocsis & Tremaine (20141. 


2.2.4 Tidal dissipation 


To investigate if tides can suppress eccentricity excitation, 
we consider the “equilibrium tide” with constant time lag to 
calculate the inner binary’s orbital evolution when the peri- 
center distance is larger than 2 Rt. Similarly to |Naoz et al.| 
(20121 and Naoz & Fabrycky (20141, we include the dif¬ 


ferential equation governing the orbital evolution following 


Eggleton et al. (1998); Egglcton fe Kiseleva-Eggleton (20011 


and Fabrycky & Tremaine (20071. For the star, we assume 


the viscous timescale is 10 yr, which corresponds to the qual¬ 
ity factor (Goldreich & Soter 19661 Q ~ 10 5 for a 10 day 
orbit (or Q ~ 4 X 10° for a 100 year orbit). 

In Figure [5] we show a representative example of the 




Figure 5. Comparison of the runs with tidal effects and the runs 
with no tidal effects. The dashed green line indicates the case 
with tidal effects and the blue lines indicates the case without 
tidal effects. The two lines are nearly identical, suggesting that 
tidal effects are negligible in these runs. The upper panel shows 
a case when a 1 Mq star orbits around a 10 7 Mq SMBH with 
ai = 0.017 pc and ej = 0.001, and is perturbed by a 10 9 Mq 
outer SMBH with <12 = l pc. The lower panel shows a case when 
a 10 Mg star orbits around a 10 7 Mg SMBH with ai = 0.035 pc 
and ei = 0.01, and is perturbed by a 10 9 Mg outer SMBH with 
<12 = l pc, e 2 = 0.7. We used the constant time lag prescription 
for the tides, and the quality factor Q was set to ~ 10 5 for a 10 
day orbit (Q ~ 4 X 10® for a 100 year orbit). 


3 SMBH-BINARY SYSTEM 

Requiring the criteria listed in §2.1[ the minimum and the 
maximum distance of the star affected by the EKL mecha¬ 
nism from the inner SMBH can be calculated. However, not 
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all stars in this region will be disrupted, since the excitation 
of the eccentricity depends sensitively on the orbital orien¬ 
tation, and the parameter region where the eccentricity can 
be excited is complicated ( Li et al.~]( )14b|. In addition, when 
the Kozai timescale is only slightly smaller than the GR 
or the NT timescale (with tx still being the smallest), the 
evolution of the inner orbit is complex. For instance, the ec¬ 
centricity of the inner orbit can be excited in configurations 
where the eccentricity cannot be excited due to the Kozai- 
Lidov mechanism alone. This excitation may be caused by 
the resonances between the NT, GR or Kozai-Lidov preces¬ 
sions (Naoz et ah 2012 1 . 

We consider the following illustrative example: mo = 
10 7 Mq, m 2 = 10 s M 0 , 02 = 0.5 pc, e 2 = 0.5. We adopt 
the isotropic stellar distribution function of equation 
assuming the stars have a solar mass, and that the ec¬ 
centricity distribution is thermal ( dN/de = 2e). We run 
large Monte-Carlo simulations, integrating the equations 
presented in Sj2] where the equations of motion for the EKL 
mechanism are given by Hamilton’s equations (eqs. (A26- 
35) in |Naoz et al. 0 13a| ), and gi = q\,ekl + gi,GR + 5i,jvt, 
32 = (/ 2 ,ekl + g 2 ,GR- We distinguish three outcomes for the 
EKL evolution: “TDE”, “scattered by the SMBH compan¬ 
ion”, and “surviving”, as explained now. 

The eccentricity of the star needs to reach very close 
to unity to cause tidal disruption. The tidal radius is Rt = 
5 x 10 -6 pc around a 10 7 Mq SMBH. We identify the TDE 
with ai(l — ei) < 3 Rt, since the stars may still be dis¬ 
rupted due to accumulated heating under the strong tide 
outside the tidal radius (Li & Loeb 2013). Since the size 
of the Hill sphere of the less massive SMBH is small (i.e. 
0.08 pc in our example), the star may reach the apocenter 
outside the Hill sphere before disruption as the eccentricity 
increases. Namely, the gravitational pull of the companion 
SMBH (m 2 ) will be larger than mo. We refer to this as a 
“scattering event” (ai(l + ei) > 02(1 — e2)(mo/(3m2)) 1/,3 ). 
Note that the secular approximation is no longer valid for 
the scattering events. Three-body integrations of the dynam¬ 
ical evolution of scattering events show that they may either 
lead to an exchange interaction, where the star is captured 
by the outer SMBH, they may cause the ejection of the star 
producing a hyper-velocity star (Samsing 2014 Guillochon 


& Loeb 20141, or they may be tidally disrupted. The scatter¬ 


ing events resulting in a capture may systematically increase 
the eccentricity distribution of stars orbiting the companion 
SMBH. For the third category, we label the stars neither 
disrupted nor scattered by the companion after 1 Gyr as 
“survivors”. 

Figure [6] shows the results of the numerical simulation 
in the final a 1 — i and a 1 — ei planes. We use open circles to 
mark stars that underwent TDEs, crosses for stars that were 
scattered by the companion, and full circles for stars that 
survived. The disruption/scattering time is color coded, and 
it indicates that most of the disruption events occur within ~ 
0.5 Myr. This corresponds to the octupole Kozai timescale, 
which is roughly 0.2 — 2 Myr for these systems at a 1 = 
0.03 — 0.08 pc. Out of all 1,000 stars between ai = 0.0275 
pc and 0.075 pc, 57 are disrupted, and 726 are scattered 
by the outer black hole. According to the stellar density 
distribution in equation (131, there are ~ 10 5 stars in this 
semi-major axis range. Normalized by the total number of 
stars in this semi-major axis range, it indicates that the tidal 


disruption rate is ~ 10 _2 /yr in the first ~ 0.5 Myr for the 
less massive black hole due to EKL, while ~ 7 x 10 4 stars 
undergo scattering events by the outer SMBH. 

Since the eccentricity of the stars with high inclinations 
are more likely to be excited, the stars with high inclinations 
are more vulnerable to tidal disruption, the final inclination 
distribution is no longer isotropic (the lower panels in Figure 
[ 7 } and the stars around the SMBH form a torus-like config¬ 
uration (see Naoz & Silk (2014]) for similar results). The 
stars with larger semi major axis have higher probability to 
be scattered when their eccentricity become excited due to 
the EKL mechanism, and thus the final distribution of stars 
surrounding the less massive black hole will be truncated 
at a larger semimajor axis. In addition, the distribution of 
the eccentricity of the surviving stars shows deviations from 
thermal distribution with a suppression of very eccentric 
stars (as expected since they get scattered by m 2 more eas¬ 
ily, and their eccentricity can be excited more easily at a 
lower inclination (Li et al. 014al). Furthermore, as shown in 
Figure [8] the stars that are closer to mo (< 0.04 pc) have 
an eccentricity distribution closer to thermal. The stars that 
are closer to m 2 (> 0.04 pc) have systematically smaller ec¬ 
centricities. The thermal distribution for closely separated 
stars (< 0.04 pc) is similar to the observed S stars in the 
center of the Milky-Way galaxy (Genzel et al. 20101, which 
shows a steeper slope. 
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Figure 6. The outcome of the evolution around a SMBH binary 
with mo = 1O 7 M0, 777-2 = 1O 8 M0, 0-2 = 0.5 pc, e 2 = 0.5. We plot 
the final i\ versus ai and ei versus ai for stars that survived, 
were disrupted, or were scattered in the simulation after 1 Gyr. 
The color code indicates the time when the star is disrupted or 
is scattered. Out of the 1,000 stars between ai = 0.0275 pc and 
0.075 pc, 57 are disrupted, and 726 are scattered by the outer 
black hole. The number of s tars in this range according to the 
distribution of equation ( | 13| > is ~ 10 5 (assuming the stars are 
1 solar mass). This suggests that the tidal disruption rate is 
10 _2 /yr in the first ~ 0.5 Myr for the less massive black hole. 


4 SMBH-IMBH SYSTEM 

Let us consider next the perturbations of a SMBH on stars 
orbiting an intermediate mass black hole (IMBH). IMBHs 
may form through runaway mergers during core collapse in 
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Figure 7. The initial distribution and the final distribution of 
the stars after 1 Gyr in our illustrative example shown in Figure[6] 
The final distribution represent the surviving stars. 



Figure 8. The final cumulative distribution of the eccentricity of 
stars in our illustrative example for mo = 10 7 Mq, m 2 = 10 8 Mq 
separated by 0.5 pc in an eccentric orbit with <:■-> = 0.5. For stars 
at distance larger than 0.04 pc, the final eccentricity distribution 
becomes shallower than that inside of 0.04 pc. 


globular clusters (Portegies Zwart & McMillan 2002 I. Since 
globular clusters sink to the galactic center through dynam¬ 
ical friction, and the disrupted globular cluster could con¬ 
tribute to most of the mass in nuclei stellar cluster for galax¬ 
ies with total mass below IO^Mq, this setup may be com¬ 
mon in the Universe(|Portegies Zwart et al.|[2006| |AntoniiiI| 


2013 Gnedin et al. 20141. Alternatively, IMBH may form 


at cosmologically early times from population III stars in 
galactic nuclei (Madau fc Rees || 2001 [) , or in accretion disks 


around SMBHs (Goodman & Tan 2004 McKernan et al. 


2012 2014). In the Milky Way center, the orbits of the S- 


stars are consistent with that caused by the dynamical in¬ 
teractions of IMBHs (Merritt et al. 20091. In addition, IRS 
13E may potentially host an IMBH, though its existence is 
controversial (|Maillard et al.|2004||Schodel et al.|2005||Fritz 


et al.||2010l. The TDE rate has been discussed by Chen & 


Liu (2013) and Mastrobuono-Battisti et al. (2014). Here, 


we consider the interactions of stars surrounding IMBHs in 
the center of galaxies with the central SMBH due to the 


hierarchical three body interactions, and consider the re¬ 
distribution of the stars as a result of the interaction. 

We set the IMBH mass to 10 4 M@ at a distance of 0.1 
pc from Sgr A* (02 = O.lpc, e 2 = 0.7, mo = 10 4 Mq and 
m 2 = 4 x 10 6 Mq). These parameters for the IMBH are al¬ 
lowed according to limits on the astrometric wobble of the 
radio image of Sgr A* (Hansen & Milosavlje vi6|[2003 ; Reid 
& Brunthaler 20041, the study of hypervelocity stars (Yu 


& Tremaine 20031, and the study of the orbits of S stars 
(G ualandris &; Merritt||2~009 ). We set the distance of stars 
to be uniformly distributed between 0.00045 pc and 0.0028 
pc. The tidal disruption radius for 1 Mg stars is 4.89 x 10~‘ 
pc. The minimum distance is set by requiring the GR pre¬ 
cession timescale to be longer than the Kozai timescale, and 
the maximum distance is set by requiring the stars to stay 
in the Hill sphere of the IMBH. Note that in this case the hi¬ 
erarchical criterion (i) in Sec. 2.1 e < 0.1, is satisfied as long 


as the stars are within the IMBH’s Hill sphere. We assume 
the distribution of the stellar eccentricity to be uniform. We 
take into account GR precession, NT precession and EKL 
at octupole order in the integration. 

In 1,000 runs, we find that ~ 40 end up in tidal dis¬ 
ruption and ~ 500 are scattered as shown in Figure [9] The 
tidal disruption/scattering time (color coded) is around 10° 
yrs. As shown in Figure |10| we predict that the surviving 
stars form a torus-like configuration (similarly to the result 


achieved by Naoz & Silk (2014) for dark matter particles). 


The predicted distribution may be resolved if the angular 
resolution of the instrument is better than that correspond¬ 
ing to the Hill sphere around the IMBH, in this case 0.07 
arcsec. This can be achieved in near infrared by the Gemini, 
VLT and Keck telescopes. In addition, the EKL mechanism 
also produces scattering events which may be responsible for 
the observed hypervelocity stars. The TDE rate may reach 
~ 10 4 / yr for a short ~ 10 5 yr duration episode after the 
globular cluster first approaches the galactic nucleus at a 
distance of 0.1 pc, assuming there are ~ 200 stars in a glob¬ 
ular cluster around an 10 4 A/q-IMBH in the EKL-dominated 


region according to the density distribution in equation (131. 
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Figure 9. The final distribution of stars surrounding a 1O 4 M0 
IMBH at a distance of 0.1 pc from Sgr A* after 100 Myr. The 
open circles represent stars that get tidally disrupted, and the 
crosses represent stars that get scattered. Both are colored ac¬ 
cording to the time of tidal disruption/scattering. We find that 
~ 50% of the stars survived tidal disruption and scattering. The 
final distribution of the star has a deficiency at high inclination 
relative to the orbital plane of IMBH. 
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Figure 10. The initial distribution and the final distribution of 
the stars after lOOMyr in our illustrative example for the IMBH, 
as shown in Figure [9] 


5 CONCLUSIONS 


SMBH binaries are natural outcomes of galaxy mergers. An 
SMBH binary may show an enhanced TDE rates due to the 
eccentric Kozai-Lidov (EKL) mechanism and chaotic three 
body interactions (Ivanov et al.|2005 Chen et al. 2009 2011 


Wegg fe Bode|2011 1 . The higher tidal disruption rates may 

in turn serve as a flag to identify closely separated black 
hole binaries on subparsec scale, which are difficult to detect 
otherwise. We focused on the effect of the EKL mechanism 
(see |Naoz et al. 2011 013al on the surrounding stars in 
SMBHB. This mechanism can excite the stars’ eccentricity 
to values very close to unity (e.g., Naoz et al. 013a|b Li| 
et al. |014a|b| ). We identified the range of physical parameters 
where EKL is important. 

We first compared the Kozai timescale with the secu¬ 
lar timescales of other mechanisms that may suppress EKL 
in galactic nuclei. These include Newtonian (NT) precession, 
general relativistic (GR) precession, resonant relaxation, two 
body relaxation, Lense-Thirring precession and orbital de¬ 
cay due to gravitational wave emission. We have found that 
for the SMBHB cases we considered, NT precession and 
GR precession may suppress EKL, especially when the inner 
SMBH is more massive than the outer SMBH (as shown in 
Figure[4|. This is consistent with the results by Naoz & Silk 


|2014| for dark matter particles around SMBH binaries, as 
well as the three body scattering experiments done by | Chen] 
et al. (20091; Wegg & Bode (20111; Chen et al. (20111, who 


observed that the tidal disruption events were dominated by 
the three body chaotic interactions rather than EKL mech¬ 
anism for stars surrounding the more massive black hole. 
However, we found that a massive outer binary allows a non- 
negligible region of parameter space where the EKL mecha¬ 
nism may operate and lead to TDEs. We also demonstrated 
that tidal effects are typically negligible for the stellar or¬ 
bital evolution (see Figure |5|. 

To illustrate the EKL effects on stars surrounding the 
less massive black hole, we ran 1,000 numerical experiments 
with different initial conditions for a star cluster surrounding 
a 10‘ Mq black hole, which is being perturbed by a 10 8 Mq 
outer black hole. We have found over ~ 50 out of the 1,000 
runs stars are disrupted in ~ 0.5 Myr. Scaled with the total 
number of stars according to equation (|13|), this corresponds 


to a TDE rate of 10 /yr for the first ~ 0.5 Myr. In con¬ 


trast, Chen et al. (20111 considered tidal disruption rates for 
stars surrounding the more massive SMBH, using numerical 
three body scattering experiments. They estimated the tidal 
disruption rate to be as high as 0.2 per year mainly due to 
three-body scattering effect^] in the first 3 x 10 5 yrs for stars 
surrounding a 10 7 Mq SMBH perturbed by an 81 times less 
massive outer SMBH. For the same SMBHB configuration, 
EKL only affects at most ~ 10 3 stars surrounding the less 
massive SMBH as shown in Figure [4] and affects at most 
~ 10 3 stars surrounding the more massive SMBH. Thus, 
EKL contributes negligibly to the total tidal disruption rate 
in this case, but EKL contributes significantly to the TDE 
rate of stars around the secondary SMBH. 

The EKL mechanism also affects the stellar distribution 
for stars surrounding the less massive SMBH. As shown in 
Figure]?] the survived stars within a particular range of radii 
are distributed in the shape of a torus ( Naoz &; Silk|]2014 |. 
In addition, a large number of stars orbiting the less massive 
black hole will be scattered by the outer black hole following 
the EKL-induced eccentricity increase. In our illustrative ex¬ 
ample, ~ 670 out of 1000 stars are eventually transferred to 
an orbit around the outer, more massive SMBH. This may 
produce hyper-velocity stars ( Guillochon fe Loeb||2014 l. 

Finally, we studied the tidal disruption of stars by an 
IMBH during mergers of globular clusters with galactic nu¬ 
clei. For an IMBH of mass 10 4 Mq at a distance of 0.1 pc 
from Sgr A*, 4% of stars get disrupted within the relevant 
distance range around the IMBH, and ~ 50% get scattered 
within 10 5 yrs. This yields a temporary tidal disruption rate 
of ~ 10 _4 /yr. Some of the scattering events may produce 
hypervelocity stars or additional TDEs. The EKL mecha¬ 
nism produces a torus-like stellar distribution for the sur¬ 
viving stars, which may be resolved by the Gemini, VLT 
and Keck telescopes in near infrared. Further investigations 
of this process using numerical scattering experiments would 
be a worthwhile in the future. 
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